##get the command line arguements
args=(commandArgs(TRUE))

#evaluate the arguments
for(i in 1:length(args)) 
	
	{
 
	eval(parse(text=args[[i]]))

	}

#Establish directories

in.dir = '/home1/99/jc152199/brt/data/'
setwd(in.dir)
out.dir = '/home1/99/jc152199/brt/data/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('gbm','SDMTools','random','sp','rgdal')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Read in data to model

model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv',header=T)

# Randomly shuffle model.data before reading into gbm()

model.data = model.data[c(sample(c(1:nrow(model.data)), nrow(model.data),replace=FALSE)),]

# Establish a list of learning rates to test

lr.list = c(.5,.1)

# Establish a for loop to create models with multiple learning rates
# Within the loop create plots to compare CV Error and Valid Error

for(lr in lr.list)
	
		{
		
		# Run gbm model
		
		brt.gbm = gbm(formula = micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, distribution = "gaussian", data = model.data, n.trees = 500, interaction.depth = 1, shrinkage = lrs, train.fraction = .1, bag.fraction = 0.5, cv.folds=10)
		
		# Determine optimal tree number using cross-validation on the model data
		
		op.tree.num = gbm.perf(brt.gbm, plot.it = FALSE, oobag.curve = FALSE, overlay = TRUE, method='cv')
		
		#Establish x and y limits for plot

		ylims = c(0,max(c(max(brt.gbm$cv.error),max(brt.gbm$valid.error))))
		xlims = c(0,500)
		
		# Reset working directory for specific learning rate
		
		setwd(paste(in.dir,lrs,'/',sep=''))
		
		# Open the .png device driver
		
		png('cverror.validerror.compare.png')

		# Set the plotting space
		
		plot(c(seq(1,500,1)), brt.gbm$valid.error, xlab = 'Iteration', ylab = 'Deviance', main = 'Error Comparison', col=1, type='n', xlim=c(xlims), ylim=c(ylims))
		
		# Create a legend for the plot
		
		legend("topright", legend=c('CV Error','Valid Error'),text.col=c(1,2),bty="n")
		
		# Plot the data
		
		lines(brt.gbm$cv.error,col=1)
		lines(brt.gbm$valid.error,col=2)

		# Flush the .png device
		
		dev.off()
		
		}
		
# End loop